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CPMC 
1. Introduction 

In these lectures we describe the constrained path Monte Carlo (CPMC) 
method for quantum many-fermion systems. We will focus on the ground- 
state CPMC algorithm[l] and its applications [2], but will also briefly de- 
scribe work in progress on a finite-temperature extension[3]. 

In previous lectures we have seen both the auxiliary-field quantum 
Monte Carlo (AFQMC)[4, 5] and Green's function Monte Carlo (GFMC) 
[6] methods. We will see here how the ground-state fermion CPMC algo- 
rithm combines the concepts of Hubbard-Stratonovich transformation and 
Slater determinants with branching random walks and importance sam- 
pling. We will demonstrate the origin of the fermion sign problem[7, 8] 
in this context and then discuss an approximate solution, the constrained 
path approximation. CPMC is free of any sign decay and its computing 
time scales algebraically with system size. The Hubbard model will be used 
as an example to illustrate the basic idea and actual implementation of the 
algorithm, although the formalism of the algorithm is more general. 

Having introduced the Hubbard model and discussed the ground-state 
CPMC method, we will then describe in some detail an application of 
CPMC. The Hubbard model has been the subject of intense theoretical ef- 
fort for the past decade, to study its ground-state properties and understand 
whether it contains the relevant electron correlations to describe high-Tc 
superconductivity [9, 10]. CPMC enabled, for the first time, ground-state 
calculations on large system sizes. On the other hand, the task of charac- 
terizing the electron pairing correlations is an extremely difficult one and 
requires great accuracy. This application illustrates both the promise and 
the challenge facing our current fermion methods. 
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In the last part we will briefly discuss some preliminary results on a 
finite-temperature (T > OK) CPMC method. The goal is to have an al- 
gorithm which enables calculations at finite temperatures free of any sign 
decay, while maintaining much of the grand-canonical formalism of the 
standard Blankenbecler, Scalapino, and Sugar (BSS) algorithm[4]. 

Why do we need another approximate fermion method? After all, we 
have fixed- node [11, 12] Green's function Monte Carlo (GFMC) (includ- 
ing diffusion Monte Carlo (DMC)[11, 13]) and projector quantum Monte 
Carlo[5, 14], for ground states, and path integral Monte Carlo[15, 16] and 
auxiliary-field quantum Monte Carlo[4] for finite-temperature properties. 
The short answer is that there are many applications that demand it. In 
fact the ground-state CPMC algorithm has already made possible various 
calculations and is seeing an increasing number of applications in different 
areas. We will also give a long answer to this question by making a rough 
table of the existing algorithms, as shown on the next page. We loosely 
categorize these fermion algorithms into continuum and lattice methods. 
The table is a synopsis to capture some basic ideas, which will provide a 
very brief review of some of the materials we have learned, and also help set 
the reference frame for the algorithms we will discuss. (Readers unfamiliar 
with AFQMC may wish to first read Section 2.1.) 

While cross-fertilization certainly has occurred, the division line be- 
tween the two classes of methods remains quite visible. Each has distinct 
strengths, which have to a large extent maintained their separation in ap- 
plication areas. Although each has seen a great deal of success, limitations 
do exist in both, as we see from the table. In this summer institute, we 
have seen considerably more "mixing" of the two, and we can hope that 
this will stimulate even more such efforts. 

2. Ground-state (T = OK) CPMC algorithm 

The ground-state constrained path Monte Carlo (CPMC) algorithm filled 
in the void (indicated by the first set of question marks in Table 1.) of a 
ground-state auxiliary-field algorithm without exponential scaling. CPMC 
is free of any "sign decay" , which is the signature of the fermion sign prob- 
lem. The required computation time therefore scales algebraically with sys- 
tem size. It is approximate, with results dependent on the trial wave func- 
tion that is used as constraint. Similar to the fixed-node approximation 
[11, 12] in configuration space, the process of CPMC solves the Schrodinger 
equation under some boundary condition defined by the trial wave function. 
But, differently from fixed-node, this boundary condition is not limited to 
real space, and can be a more "global" one, as we will discuss in detail 
below. CPMC also retains (at least to a large degree) an important advan- 
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TABLE 1. Synopsis of standard fermion quantum Monte Carlo (QMC) algorithms. 



Continuum 



Lattice 



Applications electronic structure, quantum chem- 
istry, 3 He, few-body nuclear physics 

Basis configuration space: 

R = {n, f2, . . . ,fjv}, where f; is posi- 
tion of the j-th of N fermions 

Ground-state: 



Algorithm 
Description 

P ee e- ATH 
(At > 0) 

|*(°)>: 
known 
initial 
w.f. 



Green's Function MC (GFMC)[6, 13] 
($ (o) |pp . . . p . . . -> ground state. 



■o:.-.-;. 



<)-. 



::: 



R 

Wave function amplitude is given by 
distribution of a collection of R's. Each 
application of P is realized by random 
walk for each R: branching + diffusion. 



correlated electron models, nuclear shell 
model, quantum field theory 

Slater determinant space: 

W = {01 1 <hi ■ ■ ■ j 4>n}, where each 

is a complete single-particle orbital. 



Projector QMC (PQMC)[5,14] 

(*(°)|P...O...P|*(°)> 



Compute 



(m\p...p\m) 



•(I 



)■ 



Sign problem Caused by P's that represent + and — 
amplitudes of the wave function collaps- 
ing to the same asymptotic distribution 
— a symmetric one (bosonic). 

fixed-node approximation[ll, 12] avail- ??? 
able to achieve algebraic scaling 

Observables Hard to compute (R's are orthogonal). Easy 
Finite- temperature: 



.(.-space 

P can be written as J dxp(x)B(x). 
B(x) is 1-electron operator and depends 
on auxiliary-field x. MC samples x 
("spins") to evaluate integrals. 
( | and | ) denote Slater determinants. 
( | B -» ( | ; ( | ■ | ) -> number. 

Caused by cancellations: numbers of 
paths that lead to — and + contribu- 
tions become increasingly equal as path 
grows longer (lower temperature). 



Algorithm Path-Integral MC (PIMC)[15] 

Description Feynman path integrals allow mapping 
to a classical problem of a ring-polymer 
in P-space, whose length is oc 1/T. MC 
is used to sample the "classical" parti- 
tion function (to compute observables). 



(Grand-Canonical) QMC[4] 

Similar to above. However, |$( ') is 
implicit: N is not fixed; trace over all 
possible values can be evaluated analyt- 
ically which yields expression involving 
only x. Integral over all x is again done 
by MC. (Length of path oc 1/T.) 



Sign problem restricted path-integral appr. [16] 



??? 
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tage of AFQMC compared to continuum methods, namely the ability to 
compute easily many expectation values, including off-diagonal ones such 
as electron pairing correlations. This is an important point since in many 
applications, particularly those in the study of strongly correlated models, 
it is crucial to have information on properties beyond the total energy. 

The ground-state CPMC algorithm has two main components: The first 
is to formulate the projection of the ground state as open-ended random 
walks with importance sampling, as in GFMC, although the random walks 
take place in a space of Slater determinants. The second is to constrain 
the paths of the random walks so that any Slater determinant generated 
maintains a positive overlap with a known trial wave function \ipr). Below 
we describe them in sections 2.2 and 2.3 respectively. The first is an alter- 
native way to do projector QMC which in many cases can be more efficient 
than the latter, because of the use of importance sampling. Like projector 
QMC, it is exact, but suffers from the sign problem. The second component 
is the constrained path approximation, which eliminates the sign decay, but 
introduces a systematic error in the algorithm. These two components are 
independent of each other, and can be used separately. Their combination 
is what we are calling the ground-state CPMC algorithm. 

In describing the CPMC algorithm, we will use the Hubbard model and 
often use a trivial system in our illustration. We do so to make the descrip- 
tion more pedagogical. However, it is important to think about how the 
method generalizes. This part of the notes should be used to complement 
Ref [1], in which more formal and rigorous discussions can be found. In 
addition, details of the algorithm that we will not be able to fully cover 
here can be found in those papers. 

2.1. BACKGROUND 

Here we introduce some background materials and illustrate them with a 
simple example. Some of these have been discussed by A. Muramatsu in 
his lecture[17]. 

2.1.1. The Hubbard Hamiltonian 

The one-band Hubbard model is a simple paradigm of a system of inter- 
acting electrons. Its Hamiltonian is given by 

H = K + V=-tY, (4cj<t + sW) + n ^ n ih (!) 
(ij)a i 

where t is the hopping matrix element, and c\ a and q ct are electron creation 
and destruction operators, respectively, of spin a on site i. The notation 
( ) indicates near-neighbors. The on-site Coulomb repulsion is U > 0, and 
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nia = c] a Cia is the electron number operator. We will denote the number of 
lattice sites by N, and the linear dimension by L. The numbers of electrons 
with spin a =|, j, which we denote by and N^, will be fixed in each 
calculation. That is, we will work in the canonical ensemble. We will set 
t = 1 and impose periodic boundary conditions. 

2.1.2. The Hubbard- Stratonovich transformation 

We now introduce Hirsch's discrete Hubbard-Stratonovich (HS) transfor- 
mation^]. For any positive At, which we will choose to be small, the 
following identity holds: 



where 7 is determined by cosh(7) = exp(ArC//2). For reasons that will 
become clear in the next section, we have inserted a function p{xi) in place 
of the constant factor 1/2 and will interpret p(xj) as a (discrete) probability 
density function of Xi, with Xi = ±1. 

The essence of the HS transformation is the conversion of an interact- 
ing system into many non-interacting systems living in fluctuating external 
auxiliary-fields, and the summation over all such auxiliary-field configura- 
tions recovers the many-body interactions. More specifically, the goal is to 
write the many-body operator e~ ArH in terms of one-electron operators. In 
Equation (2), the exponent on the left, which comes from the interaction 
term V, is quadratic in n, indicating the interaction of two electrons. The 
exponents on the right, on the other hand, are linear in n, indicating two 
non- interacting electrons in a (common) external field characterized by X{. 

The HS transformation above is one special case for the Hubbard in- 
teraction. Since our purpose is to illustrate the CPMC algorithm with an 
example, we will not dwell on more general forms of HS transformations [19] 
here, except to state that they exist for various other forms of interactions. 
In general, after an HS transformation is applied, the propagator e~ ArH 
can be written as 



where p(x) is a probability density function and B(x) is a one-electron 
operator whose matrix elements are functions of the many-dimensional 
auxiliary-fields x. 

2.1.3. A specific system to illustrate things 

We now look in some details into the ground-state of a simple 2x2 system 
of the Hubbard model to help explain the "language" we will use when 
describing CPMC later. We will consider iVf = 2 and iVj_ = 1. We label the 
four sites 1 thru 4 such that site 1 is near-neighbors with sites 2 and 3. 



e -Ar[/n iT n a _ 



-ArC"(n iT -t-nii)/2 / "UTS. («»t ) 



(2) 





(3) 
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First let us examine the trivial case of free electrons, i.e., we set U = 0. 
We can write down the 1- electron Hamiltonian matrix, which is of dimen- 
sion 4x4: 



H = 



( 





-1 


-1 





\ 




-1 








-1 




-1 








-1 




\ 





-1 


-1 








The eigenstates of H can be obtained by direct diagonalization. With these 
eigenstates, we immediately obtain the ground-state wave function |^o) of 
the 3-electron system from the Pauli exclusion principle: 



( 0.5 -0.58 \ 

0.5 0.40 

0.5 -0.40 

V 0.5 0.58 J 



/ 0.5 \ 

0.5 

0.5 
V 0.5 / 



where the first matrix contains two single-particle orbitals (two columns) 
for the two f electrons and the second matrix contains one single-electron 
orbital for the one J. electron. Each single-electron orbital is an eigenvector 
of H. Note the second and third lowest single-electron states are degenerate, 
causing the 3-electron system to be open shell and \ipo) to be degenerate. 
We have simply picked one particular linear combination for the second 
f-electron in |^o) above. 

An object of the form of |V>o) is of course nothing more than a Slater 
determinant. For example, the amplitude of the configuration \R) = | j 
It), i.e., two | electrons on sites 3 and 4 and the one j electron on site 1, 
is given by 

Wo) = det ( ° Q j ^ ) • det ( 0.5 ) . 
That is, more formally, 



\R) = 



1 





\ 













1 







® 


V o 


1 


) 


Vo; 



and (R\ipo) = det(i?^ • ^o)) where R and denote the matrices corre- 
sponding to \R) and |^o)j respectively. In general, the overlap of two Slater 
determinants 

(<j)\<t>') = det($t . $') (4) 

is a number. 
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For this non-interacting system, an alternative (albeit indirect and in- 
deed circular) way of obtaining \ipo) is by the power method. From the 
eigenvalues and eigenvectors of K, we can easily construct the matrix for 
e ~ArH ^ w \ l [ c \ 1 tn e same structure as H above (i.e., a 4 x 4 matrix). 
Denote this matrix by Bk- With an arbitrarily chosen initial Slater de- 
terminant |?// )) (with non-zero overlap with |V>o))) we can then repeatedly 
apply e~ ArH , which means multiplying both the 4 x 2 j matrix and the 4x1 
I matrix by Bk- As indicated in Table 1, this leads to \ipo) asymptotically, 
i.e., 

|V (n+1) ) oce- ArH |V {n) ) (5) 

gives IV'o) as n ^ oo. 

Now suppose we turn on the interaction U. The first approach of di- 
rectly diagonalizing H is the method of exact diagonalization, which scales 
exponentially. The power method of Eq. (5), on the other hand, can still 
apply if we can write e~ ArH in some one-electron form. The HS trans- 
formation does just that. Assuming At is small and applying the Trotter 
break-up, we have 



-AtH 



( e 



7x1 













































\ 



Bk 



J 



V 



-7m 







-■yx 2 

















-7x4 



Bk, 



where x = {x±, X2, X3, X4}. This is just Eq. (3). Note that B(x) has an 
\ and a J. component, each of which is a 4 x 4 matrix. Applying each 
B(x) to a Slater determinant means precisely the same as in the non- 
interacting case (with Bk <8> Bk)- In other words, B(x) operating on any 
Slater determinant \<f>) simply involves matrix multiplications for the | and 
j components separately, leading to another Slater determinant \<j)'): 



b') = B{x)\ct>). 



(6) 



We now have all the basics to understand the projector QMC algorithm 
described in Table 1, and to move on to discuss CPMC. 



2.2. RANDOM WALK IN SLATER DETERMINANT SPACE 

The first component of the CPMC algorithm is the reformulation of the 
projection process as branching, open-ended random walks (a la GFMC) 
in Slater determinant space. 
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We start from the projection process in Eq. (5). Using (3), we write it 

|V> (n+1) > = J dxp{x)B{x)\^ n) ). (7) 

In the Monte Carlo realization of this iteration, we represent the wave 
function at each stage by a finite ensemble of Slater determinants, i.e., 

|Y> (n) >oc£l^ n) >> (8) 
k 

where k labels the Slater determinants and an overall normalization factor 
of the wave function has been omitted. Analogous to GFMC, the Slater de- 
terminants will be referred to as random walkers. We note that, contrary to 
the walkers in GFMC, the walkers here are non-orthogonal. Further, their 
characterization is not as simple as that of a point, because various combi- 
nations of single-particle orbitals can lead to the same Slater determinant. 

The iteration in (7) is achieved stochastically by Monte Carlo (MC) 
sampling of x: 

J dxp(x)B(x)\^); (9) 

that is, for each random walker we choose an auxiliary-field configuration 
x from the probability density function p(x) and propagate the walker to a 

new one via |</4™ +1 ^) = B(x)\cj)^). We repeat this procedure for all walk- 
ers in the population. These operations accomplish one step of the ran- 
dom walk. The new population represents \ip( n+1 ^) in the sense of (8), i.e., 
|^(n+i)^ ^ \<j)^' +l ' > ). These steps are iterated indefinitely. After an equi- 
libration phase, all walkers thereon are MC samples of the ground-state 
wave function |^o) and ground-state properties can be computed. 

The algorithm we have described so far, while correct, is not efficient, 
because the sampling of x is completely random with no regard to the po- 
tential contribution to the ground-state wave function. In order to improve 
the efficiency of (7) and make it a practical algorithm, an importance sam- 
pling^] scheme is required. To further motivate the need for importance 
sampling and develop some insights on how to proceed, we consider the 
so-called mixed estimate of the ground-state energy. As in GFMC, this 
is given by Eq = {i{jt\H\iPo)/(iPt\iPo}- Estimating Eq requires estimating 
the denominator by S<a(V"t in which \<p) denotes random walkers after 
equilibration. Hence, the term (ipxl^) plays the role of the Boltzmann distri- 
bution in statistical physics, while the denominator resembles the partition 
function. Our current way of sampling corresponds to sampling completely 
randomly configurations from the ensemble, which would clearly lead to 
large statistical fluctuations in the evaluation of the partition function. To 
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improve the situation, it is natural to draw the analogy with the standard 
practice of Monte Carlo in statistical physics and try to sample \<j>) ac- 
cording to (V>t !</>)• This way, each walker would ideally contribute equally 
to the estimate, and Eq would be given by the average of the quantity 
(iPt\H\(/)) / (iPt\</>) w hh respect to the sampled walkers. For a good \ipr), 
fluctuations in this quantity are small and we expect the statistical error 
to be much reduced. 

With importance sampling, we iterate a modified equation with a mod- 
ified wave function, without changing the underlying eigenvalue problem of 
(7). Specifically, for each Slater determinant we define an importance 
function 

o T {4>) = ^ T \<p), (io) 

which estimates its overlap with the ground-state wave function. We can 
then rewrite equation (7) as 

|Vi (n+1) ) = J dxp(x)B(x)\i>^), (11) 

where the modified "probability density function" is 

p(x) = , , J p{x). (12) 

As expected, p(x) is a function of both the current and future positions in 
Slater-determinant space. Further, it modifies p(x) such that the probabil- 
ity is increased when x leads to a determinant with larger overlap and is 
decreased otherwise. It is trivially verified that equations (7) and (11) are 
identical. 

In the random walk, the ensemble of walkers { \4>^) } now represents 
the modified wave function: \4>^) oc J2k l^)- The true 

wave function is 

then given formally by 

I^VEI^VOt^), (13) 



although in actual measurements it is |i/A n )) that is needed and division by 
Ot does not appear. The iterative relation for each walker is again given by 
(9), but with p(x) replaced by p(x). The latter is in general not a normalized 
probability density function, and we denote the normalization constant for 
walker k by N(<j>^) and rewrite (9) as 

l^ n+1) ) - N(^) J dx-^B(x)\^). (14) 
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This iteration now forms the basis of the CPMC algorithm. As in the DMC 
method, it is convenient to associate a weight w k with each walker, which 
can be initialized to unity. With the introduction of the weig ht, 

Sfc^fc^l^jfc 1 ^)- F° r eacn walker \4>^}, one step of the algorithm is then as 
follows: 

1. sample a x from the probability density function p(x) / 'N {4>^' ') , 

2. propagate the walker by B(x) to generate a new walker, 

3. compute a weight u/[ n+1 ^ = N((f>^) for the new walker. 

To better see the effect of importance sampling, we observe that if 
IV't) = IV'o)) the normalization J p(x)dx is constant. Therefore the weights 
of walkers remain a constant and the random walk has no fluctuation. Fur- 
thermore, we refer again to the estimator for E$. With importance sam- 
pling, the denominator becomes the sum of weights w, while the numerator 
is Y J <f ) {i'T\H\4>)w r f > / (i^t^), where again \4>) denotes walkers after equilibra- 
tion. As \iPt) approaches \ipo), all walkers contribute with equal weights to 
the estimator and the variance approaches zero. 

We stress again that the importance sampling transformation does not 
alter the underlying iterative equation in (7). The results of the calculations 
are not affected; different choices of the importance function only affects 
the efficiency of the algorithm and thus the statistical error. 

2.3. THE SIGN PROBLEM AND THE CONSTRAINED PATH 
APPROXIMATION 

The second component of the CPMC algorithm is the constrained path 
approximation to eliminate the sign decay. 

As a side- product, the reformulation of the projection process described 
in 2.2. provides a clear picture of the origin of the sign problem. Below we 
first elaborate on this picture, from which the constrained path approxima- 
tion emerges naturally. 

The sign problem occurs because of the fundamental symmetry between 
the fermion ground-state |^o) an d its negative — |^o) [20, 21]. For any ensem- 
ble of Slater determinants {\<fi)} which gives a Monte Carlo representation 
of the ground-state wave function, this symmetry implies that there exists 
another ensemble {— \<p)} which is also a correct representation. In other 
words, the Slater determinant space can be divided into two degenerate 
halves (+ and — ) whose bounding surface TV is defined by (ipo\(fr) = 0. This 
surface is in general unknown. 

In some special cases, such as the particle-hole symmetric, half-filled 
one-band Hubbard model, symmetry prohibits any crossing of M in the 
random walk. The calculation is then free of the sign problem [18]. In more 
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general cases, however, walkers do cross M in their propagation by e~ . 
The sign problem then invariably occurs. 

In Fig. 1, we illustrate the space of Slater determinants by a one- 
dimensional (horizontal) line. The node M is the (red) dot in the middle. 
Imaginary time (or n) is in the vertical direction and increases as the arrow 
suggests. That is, as the walker moves in the horizontal line, we stretch 
out continuously "snapshots" of its position along the vertical direction. 
Now we follow an initial Slater determinant. With no loss of generality, we 
assume it has a positive overlap with \ipo)- At time n = it is indicated by 
the (green) dot on the right. As the random walk evolves, the walker can 
reach the node, which is the (red) vertical line. At the instant it lands on 
AT, the walker will make no further contribution to the representation of 
the ground state, since 

(V>o|0>=O => (Vo|e- rH |0)=Oforanyr. (15) 

Paths that result from such a walker have equal probability of being in 
either half of the Slater determinant space. A few of these possible paths 
are shown by dashed lines. Computed analytically, they would cancel and 
not make any contribution in the ground-state wave function, as indicated 
by their symmetric placement with respect to the node line. But since the 
random walk has no knowledge of jV, these paths continue to be sampled 
(randomly) in the random walk and become Monte Carlo noise. Only paths 
that are completely confined to the right-hand side, as shown by the solid 
(green) line, will lead to contributions to the ground state, but the relative 
number of such confined paths decreases exponentially with n. Asymptoti- 
cally in n, the Monte Carlo representation of the ground-state wave function 
consists of an equal mixture of the + and — walkers, regardless of where the 
random walks originated. The Monte Carlo signal is therefore lost. The de- 
cay of the signal-to-noise ratio, i.e. the decay of the average sign of (V"r|</>), 
occurs at an exponential rate with imaginary time. 

To eliminate the decay of the signal-to-noise ratio, we impose the con- 
strained-path approximation. Fahy and Hamann first used [20] such a con- 
straint in the framework of the standard AFQMC method. However, the 
non-local nature of such a constraint proved difficult to implement effi- 
ciently in the "path-integral-like" projector QMC scheme. Here, with the 
random walk formalism, the constraint only needs to be imposed one time- 
step at a time and is extremely simple to implement. It requires that each 
random walker at each step have a positive overlap with the trial wave 
function \ipr): 

(VtI^Vo. (16) 

This yields an approximate solution to the ground-state wave function, 
= J2(f> \4>)i in which all Slater determinants \<f>) satisfy (16). 
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Figure 1. Schematic illustration of the sign problem. The top line represents Slater 
determinant space; the dot represents the "node" TV, where a determinant is orthogonal 
to the ground state | Vo) - As the projection continues (increasing n), Slater determinants 
undergo random walks, tracing out "paths" as shown. When a walker reaches M, its 
future paths will collectively cancel in their contribution to |V>o), indicated by the sym- 
metric distribution of dashed paths about the nodal line. The Monte Carlo sampling, 
with no knowledge of this cancellation, continues to sample such paths randomly. The 
relative number of paths with real contributions (solid paths) decreases exponentially as 
n increases. 

From (15), it follows that the constrained path approximation becomes 
exact for an exact trial wave function \iPt) = IV'o)- The overall normaliza- 
tion of walkers remains a constant on the average, since the loss of walkers 
at Af is compensated by the branching of walkers elsewhere; that is, the 
eigenvalue problem with Af as boundary has a stable solution. 

To implement the constrained-path approximation in the random walk, 
we redefine the importance function by (10): 



This prevents walkers from crossing the trial nodal surface Af and entering 
the "— " half-space as defined by {tpr}- In the limit At — > 0, (17) ensures 
that the walker distribution vanishes smoothly at Af and the constrained- 
path approximation is properly imposed. With a finite At, however, p(x) 
has a discontinuity at Af and the distribution does not vanish. We have 




(17) 
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found this effect to be very small for reasonably small imaginary-time steps 
At . Nonetheless, we correct for it by modifying p(x) near N so that it 
approaches zero smoothly at N. 



2.4. ADDITIONAL ALGORITHMIC ISSUES 
2.4.1. Computing expectation values 

Given the ground-state wave function from CPMC, where the super- 
script V indicates that the wave function is obtained under the constrained- 
path approximation, we can compute the ground-state energy with the 
mixed estimate [6]: 

pmixed _ (jHgMi) (m x 

EcPMC ~ ~wm~ m (18) 

This can be implemented in a way similar to GFMC. However, contrary to 
fixed-node in GFMC, the ground-state energy computed from the mixed 
estimate is not an upper bound[22]. (The argument given in Ref. [1] on 
variational properties is incorrect.) The origin for the disappearance of the 
upper bound property in the mixed-estimate is the non-orthogonal nature 
of the Slater determinant space. In contrast with configuration space in 
fixed-node GFMC, where points on the node have no overlap with either 
\4>t) or the fixed-node solution, a point in the Slater determinant space 
which has no overlap with \i/jt) is not necessarily orthogonal with |V>o)- 
Eq. (18) is not equivalent to the estimator (V>o|#|V , o)/(V'ol^ ; o}> which is 
variational. 

In the Hubbard model, the effect of this appears to be small[22]. In 
Table I, we show several cases we could find where the computed ground- 
state energies by CPMC fell below the exact results. The system is 4 x 4, 
with 5 | and 5 J. electrons, which corresponds to a special case of a closed- 
shell system in the non-interacting limit. The free-electron wave function 
was used as the trial \iPt)- As we see, the amount by which the mixed 
estimates fall below the exact results is no more than 0.1% across a wide 
range of interactions. For U = 8, an "almost" free-electron wave function 
was used as \4>t) m Ref. [1]> namely an unrestricted Hartree-Fock wave 
function generated with a small U. That \ipr) led to a ground-state energy 
above the true valuefl]. 

Several possibilities exist to correct for the non-variational nature of the 
mixed estimate, at least in principle, and make it an upper bound[22]. For 
example, we could compute {iPq\H\tPq) / (V'o IV'o) directly with a method sim- 
ilar to forward walking[23, 24] in GFMC. That is, we can create a separate 
walk for the left-hand wave function (ipfi | , which propagates from (ipr I with 
the constraint properly implemented, and then compute the overlap with 
populations in the regular (right-hand) walk. It would clearly be inefficient 
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TABLE 2. Examples showing the lack of upper 
bound property in the mixed- estimate of CPMC. 
The system is 4 x 4, with 5 | and 5 I electrons. 
Calculations were done for three values of the in- 
teraction strength U, 4, 8, and 20. The CPMC re- 
sults are obtained from the mixed estimate; statis- 
tical errors are in the last digit and are indicated 
in parentheses. These are compared with results 
from exact diagonalization. 



u 


4 


8 


20 


CPMC 


-19.582(5) 


-17.519(2) 


-15.460(3) 


exact 


-19.580 


-17.510 


-15.452 



to randomly sample the left-hand walks, whose overlap with the regular 
walks would fluctuate greatly. We have devised an "importance sampling" 
scheme that both imposes the constraint with \iPt) and attempts to in- 
corporate the necessary correlations between the left-hand random walker 
and its matching right-hand walker. We take an importance function 

Ofp{cf)') = \{<fi'\4>^)\ + a((j)'\ipT}, where a is a parameter and \4>u^) is the 
walker with which ((/>' | will be matched in the regular population represent- 
ing the right-hand wave function. While this offered a significant improve- 
ment over naive sampling, we found it still prone to large fluctuations. With 
the variational correction term to the mixed estimate being very small, the 
statistical accuracy obtained by this scheme did not appear adequate to 
be of practical use for large systems. We note that, although this way of 
forward walking is somewhat similar to the back-propagation technique we 
discuss below, they are not equivalent. The latter tends to be much more 
accurate statistically. However, as we see below, it does not maintain the 
correct sense of direction in imposing the constraint and does not yield 
strictly {i[)q\. 

The back-propagation (BP) technique enables the computation of the 
expectation value of an operator O that does not commute with H. As we 
have mentioned, the essence of BP also comes from the forward walking 
technique in the GFMC method: 

(0)Bp=Km (^Et^JM. (19) 

r^oo (^ T exp(-Tii c )|V>o) 

A subtle distinction, however, exists between back-propagation and for- 
ward walking. In back-propagation, (V>Texp(— tH c )\ = (ipr\ exp(— tH c ) is 
restricted to "constrained" paths, i.e., those paths that do not violate the 
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constraint in the original forward direction exp(— At H c )\?Pq). In fixed-node 
GFMC a path in configuration space has no sense of direction with respect 
to the nodal surface. In the CPMC method, however, there is a sense of 
direction: a set of determinants along the path of a random walk which 
does not violate the constraint at any step when going from right to left 
may indeed violate it any even number of times when going from left to 
right. Because of this sense of direction, expression (19) may not yield the 
true expectation value (V ; ol^lV ; o)/(V , ol , /'o)- However, since \ipo) is itself ap- 
proximate, this issue is not crucial. (O)bp does have the correct limiting 
behavior, remaining exact for an exact trial wave function. 

It is relatively simple to implement the back-propagation scheme on 
top of a regular CPMC calculation. The key is that once a path has been 
sampled, the single-particle propagator can be easily constructed. Propa- 
gating a Slater determinant with it, either forward or backward, simply 
means matrix multiplications. We choose an iteration n and store the en- 
tire population { \<p^} }• As the random walk proceeds, we keep track of 
the following two items for each new walker: (1) the sampled auxiliary-field 
variables that led to the new walker from its parent walker and (2) an in- 
teger that labels the parent. After an additional m iterations, we carry out 
the back-propagation: For each walker I in the (n + m) th (current) popula- 
tion, we initiate a determinant (ipxl and act on it with the corresponding 
propagators, but taken in reverse order. The m successive propagators are 
constructed from the stored items, with exp(— AtK/2) inserted where nec- 
essary. The resulting determinants (<^ m ^ | are combined with its parent from 
iteration n to compute (0)bp, in a way similar to the mixed estimator (18). 
The weights are given correctly by ^ n+m * ) ; due to importance sampling in 
the regular walk. Starting from another iteration n', this process can be 
repeated and the results accumulated. 

We note again that, because the random walkers in CPMC are full 
Slater determinants, the algorithm lends itself well to the calculation of 
expectation values. The non-orthogonal nature of the determinant space 
ensures that overlaps between two walkers can be computed. This is to be 
contrasted with "continuum" methods in Table 1, where off-diagonal expec- 
tations are in general difficult to compute, because it is difficult to sample 
two group of random walkers whose overlap is well behaved statistically. 

2.4.2. population control and bias correction 

In CPMC walkers carry weights or, equivalently, they branch. The proce- 
dure we use to control the population is similar to that used in many GFMC 
calculations. First, a branching (or birth/death) scheme is applied, in which 
walkers with large weights are replicated and ones with small weights are 
eliminated by probabilities defined by the weights. There exist various ways 
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to do this[25, 24, 26], with the guideline being that the process should not 
change the distribution statistically. In general, how this step is done only 
affects the efficiency of the algorithm, but does not introduce any bias. 

Branching allows the total number of walkers to fluctuate and possibly 
become too large or too small. Thus as a second step, the population size 
is adjusted, if necessary, by rescaling the weights with an overall factor. 
Re-adjusting the population size, i.e., changing the overall normalization 
of the population, does introduce a bias[25, 24]. We correct for this bias by 
carrying the m most recent overall rescaling factors and including them in 
the estimators when computing expectation values. In the calculation we 
keep a stack which stores the m latest factors, (j = l,m). Suppose 
that at the current step the total number of walkers exceeds the pre-set 
upper bound. We modify the weight of each walker by a constant factor 
/ < 1 which reduces the population size to near the expected number. 
We then replace the oldest in the stack by /. Whenever we compute 
expectation values, we multiply the weight of each walker by 1/ f\j In 
our calculations on the Hubbard model, m is typically between 5 and 10. 
As we include more such factors, i.e., increasing m, the bias is reduced, 
but the statistical error increases. On the other hand, as we reduce m, the 
statistical error becomes smaller, but the bias increases. The choice of m 
is thus a compromise between these two. Another common approach to 
eliminate or reduce bias, which is in general less efficient, is to do several 
calculations with different (average) population sizes and extrapolate to the 
infinite population limit. 

Clearly, the schemes we have described to correct for bias have some 
arbitrariness and further improvements are possible. But this should not 
be a major factor in the calculation. In cases where the bias is substantially 
larger than can be handled by correction schemes in this spirit, it is likely 
more productive to attempt to improve the importance sampling and the 
algorithm, rather than details of the bias correction scheme. 

2.4.3. re-orthogonalization a Slater determinant 

Repeated multiplications of B(x) to a Slater determinant, either in the 
regular walk or in BP, lead to a numerical instability, such that round-off 
errors dominate and \4>^) represents an unfaithful propagation of |<^}- 
This instability is well-known in the AFQMC method and is controlled [27, 
28] by a numerical stabilization technique that requires the periodic re- 

orthonormalization of the single particle orbitals in \4>^)- 

We use the modified Gram-Schmidt procedure[28, 27] to stabilize the 
Slater determinants. The procedure is further simplified because of the ran- 
dom walk formalism in CPMC. For each walker \<j>), we factor the matrix 
<3? as $ = QR, where Q is a matrix whose columns are a set of orthonor- 
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mal vectors representing the re-orthogonalized single-particle orbitals. R is 
a triangular matrix. Note that Q contains all the information about the 
walker \<f>). R, on the other hand, only contributes to the overlap of \4>). 
Since the overlap is already represented by the walker weight, R does not 
need to be carried along and can be discarded. 

It is interesting to note that this instability originates from the tendency 
for each walker (Slater determinant) to collapse to a bosonic ground state. 
In GFMC this actually happens and is the first cause for the sign problem. 
Here the instability can be eliminated by the re-orthonormalization proce- 
dure because a walker is explicitly a Slater determinant, not a single point. 
The procedure to stabilize the Slater determinants is like analytically can- 
celing walkers in GFMC [21], which requires a large density of walkers and 
does not scale well with system size. CPMC, like AFQMC, automatically 
imposes the antisymmetry. This would suggest that the sign problem is re- 
duced in this formalism compared to GFMC. One trivial observation which 
is consistent with this is the case of a non-interacting system, where CPMC 
does not suffer from the sign problem and has zero- variance, while GFMC 
does have a sign problem and still requires the fixed-node approximation. 

2.5. ILLUSTRATIVE BENCHMARK RESULTS 

Here we show a few benchmark results to illustrate characteristics of the 
CPMC algorithm. These results are for the two-dimensional Hubbard model. 
Ref. [1] contains more results and a more extensive account of the bench- 
mark calculations that have been performed. 

In Figure 2 we show the computed Eq/N at U = 4 for various lat- 
tice sizes and electron filling (n) = (iVf + N±)/N. Also shown is available 
AFQMC data [29]. We see that the CPMC results are in good agreement 
with the AFQMC data. For example, in the case of 8 x 8 lattice with 25 | 
and 25 j electrons, the CPMC result lies approximately 0.4% ± 0.2% above 
the AFQMC result. The behaviors of the statistical errors in these two al- 
gorithms are of course completely different. In AFQMC the sign problem 
causes the variance to increase exponentially with N and the projection time 
nAr, while the CPMC method is free of any sign decay and exhibits power 
law scaling with N. For large systems (10 x 10 and beyond), the CPMC 
error bars are 30-50 times smaller than those of AFQMC, as indicated for 
the 12 x 12 system. 

A more stringent test of the algorithm is the calculation of correlation 
functions. In Table 3, we show computed expectation values for the system 
of 4 x 4 with 7 f 7 [ electrons, which is an open-shell case with a severe 
sign problem. The one-body density matrix is the expectation value of the 
Green's function element: p(\) = (cqC\), where the vector 1 = (l x , l y ) denotes 
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E/N 




-1.25 

0.50 0.60 0.70 0.80 0.90 1.00 

<n> 

Figure 2. Computed energies per site vs. electron fillings from CPMC, compared with 
available AFQMC and exact diagonalization data. The lines are to aid the eye. Curves 
for L > 4 are shifted up. Half-filling is (n) = 1. For 4x4, the solid curve is from exact 
diagonalization [30]. At (n) = 0.875 (7 | 7 |), the CPMC (with error bar in the last 
digit) and exact numbers are shown, together with the variational energy from the trial 
wave function {ipr) in CPMC. Numbers are also shown for the AFQMC (lone point) and 
CPMC values for 12 x 12 with 61 | 61 |. Note that the CPMC error bar is smaller by a 
factor of 30 in this case. 

the position of a site (vs. i, which we have used so far and which is the 
site label). The spin density structure factor is 

S(k x ,k y ) = S(k) = l/JVX;exp(tk-l)(soBi), (20) 

l 
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TABLE 3. Computed expectation values and correlation functions 
from CPMC for a 4 x 4 lattice with 7 | 7 j electrons and (7 = 4, 
compared with exact results. Results are shown for two different trial 
wave functions \ipTi) and \ipT2)- Both are unrestricted Hartree-Fock 
wave functions, but |i/ti) was obtained with a U of 0.1, while |?/>t2) 
with U = 4. Exact diagonalization results are from Ref. [30]; num- 
bers in parentheses indicate either the range of values due to the 
ground-state degeneracy or uncertainties in extracting numbers from 
a graph. Statistical errors on the last digit of the CPMC results are in 
parentheses. E k is the kinetic energy, p(l x ,l y ) the one-body density 
matrix, S the spin density structure factor, and n k the momentum 
distribution. 







E k 


P(2,2) 




n fc (7r/2,0) 


variational 


\lpTl) 


-24.0 


1.654 


0.625 


1.0 




\lpT2) 


-21.88 


-0.0602 


4.39 


0.941 


CPMC 


\lpTl) 


-21.44(2) 


-0.051(1) 


2.90(1) 


0.92(1) 




\lpT2) 


-21.39(8) 


-0.049(1) 


2.92(2) 


0.92(1) 


exact 




-21.39(1) 


-0.051 


2.16(2) 


0.93(1) 



where si = — is the spin at site 1. We show results from CPMC 
simulations with two different trial wave functions. Both are unrestricted 
Hartree-Fock wave functions, but was obtained with a U of 0.1, while 

\1pT2) with [7 = 4. The calculation with \ipTi) has much less fluctuation, 
even though \1pT2) has a lower variational energy 1 . We see that the two trial 
wave functions yield very different variational estimates, but their CPMC 
results are consistent and in good agreement with exact results. For ex- 
ample, in the free-electron-like \ipri), the momentum distribution is a step 
function, and the k = (1, 0) state is completely occupied (n& = 1), but even 
with this trial wave function the CPMC method still gives the correct occu- 
pation of 0.92(1). The spin structure factor S(ir, ir) is the only case where 
CPMC does not produce the exact result. However, even with the free- 
electron-like IV'ti)) which severely underestimates the anti- ferromagnetic 
correlation, CPMC correctly recovers the physics and predicts the presence 
of a strong peak of S'(k) at (tt, it). 

The fact that CPMC shows consistency with different choices of \i/jt) is 
reassuring. Such consistency checks are crucial in real applications where 
benchmark data is scarce. The fact that CPMC appears to be able to 

1 This trend appears to be rather general: free-electron-like wave functions tend to be 
better importance functions than unrestricted Hartree-Fock wave functions. This likely 
reflects to what extent translational invariance is maintained. 



20 



SHIWEI ZHANG 



accurately calculate correlation functions in these systems is significant. 
Correlation functions are often much more difficult to compute than the 
energy. It is here that the difference in scaling behaviors between CPMC 
and AFQMC becomes much more magnified. The AFQMC algorithm tends 
to break down for system sizes much smaller than indicated in Fig. 2. The 
fixed-node GFMC also encounters difficulties in computing certain correla- 
tion functions because forward walking does not yield proper overlaps, as 
we discussed earlier. 

In Table 4, we show a comparison of the computed ground-state energies 
from CPMC and from lattice fixed-node GFMC [31]. We see that, in addi- 
tion to the ability to compute correlation functions conveniently, CPMC 
does better even for the energy. (The particular system is an easy case for 
CPMC and the results may not be typical of the Hubbard model.) At small 
U, CPMC is expected to perform better because the linear combination of 
mean-field states provides an effective description of the system. On the 
other hand, at very large U, we expect that, at least in the current way of 
Hubbard-Stratonovich transformation, the statistical error in CPMC will 
be larger than GFMC. More complete studies are necessary to characterize 
where the "cross-over" occurs. 



TABLE 4. Comparison of CPMC with lattice fixed-node (FN) 
GFMC for 4x4, with 5 | and 5 J. electrons. Shown are computed 
ground-state energies per site. The CPMC result was obtained 
with a free-electron \iPt), as was the first number for fixed-node 
GFMC. The FN GFMC result "Gutzwiller" was obtained with a 
better \tpr) which included a Gutzwiller factor. Statistical errors 
are in the last digit and are indicated in parentheses. Fixed-node 
results are from Ref. [32]. 

FN GFMC FN GFMC (Gutzwiller) CPMC exact 
-f. 2186(4) -1.2201(4) -1.2239(3) -1.2238 



3. Application: Is the Hubbard model the right one for high-T c ? 

Here we describe an application of the CPMC algorithm to study the 
ground state of the one-band Hubbard Hamiltonian of Eq. (1). In par- 
ticular, we will focus on the electron pairing correlation in the d x 2_ y 2 -wave 
channel. Most of the data in this section is from Ref. [2]. 

The Hubbard model is a simple microscopic description of a many- 
electron system which includes both band-structure and electron interac- 
tions. Since the discovery of high-T c superconductivity, it has been widely 
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believed that this model contains essential features relevant to the prop- 
erties of the Cu02 planes in the cuprate superconductors. Viewed as the 
basis for further progress on understanding the mechanism for high-T c , the 
model has been the subject of intense theoretical activity [9]. Various calcu- 
lations predict [10] that the doped model exhibits an attractive interaction 
between pairs; the extended s- and d x i_ y i- symmetries of this attraction 
are consistent with the likely symmetries of the experimentally measured 
superconducting gap [10]. Yet unobserved, however, is convincing evidence 
that the attractive interaction leads to a ground state with off-diagonal 
long-range order [9, 10]. In other words, we do not have the complete an- 
swer to whether the ground state of this model has a condensate of electron 
pairs and, if yes, what their nature is. 

Numerical approaches have been a promising tool for answering some 
of the questions on the Hubbard model. The AFQMC method in particular 
has been applied extensively and has seen a great deal of success [33, 9]. 
However, due to the sign problem, it has typically been limited to relatively 
small system sizes, high temperatures, and selected electron fillings. CPMC 
largely eliminates these difficulties, although we must bear in mind the 
approximate nature of the constrained path approximation. 

The pairing correlation function we computed is defined as 

D a {\) = (At (1)A Q (0)), (21) 

where a indicates the nature of pairing. The pair-field operator at site 1 is 
Aa(l) = E«5/a(^)[ c lT c l+<5|- c U c l+<5T]) where 5 is (±1,0) and (0, ±1). For the 
extended s-wave (s*-), f s *(S) = 1. For d-wave, fd(S) is 1 when 5 = (±1,0) 
and —1 otherwise. 

Because the algorithm has an uncontrolled approximation and because 
the electron pairing correlations are extremely difficult quantities to char- 
acterize, we have carried out benchmark calculations [2, 1] and many self- 
consistency checks. Due to the limited availability of data for benchmark, 
the latter becomes even more important. In Fig. 3, we show another such 
example, for a fairly difficult case of 8 x 8 with 27 f 27 j electrons. 

Fig. 4 shows the d x 2_ ?; 2-wave pairing correlations for 12 x 12 and 16 x 16 
lattices at U = 2 and 4. The electron filling is 0.85, which corresponds to 
closed-shell cases, with = N t = 61 for 12 x 12 and jV T = = 109 
for 16 x 16. The free-electron wave function was used as \^t) m CPMC. 
Also shown in each graph is the D^(l) for the non-interacting system. We 
see that the computed d x 2_ y 2-wave pairing correlations at both values of U 
are bounded from above by the U = result. Further, consistent with this 
observation, we see that at large distances D ^(1) decreases as U increases 
from 2 to 4. In the 16 x 16 case, the long-range correlations vanish as the 
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Figure 3. Sensitivity of CPMC results to the choice of trial wave function \^t)- Shown 
is the near-neighbor d-wave pairing correlation as a function of pair separation |1|. The 
short distance portion is omitted. This system corresponds to a difficult open-shell case. 
The two trial wave functions are both unrestricted Hartree-Fock wave functions. The first, 
\ipTi), was obtained with a very small interaction strength U, resembling a free-electron 
wave function, while \ipT2) was generated with the true U of 4. The CPMC pairng cor- 
relations, with \ipTi) and \1pT2) as trial wave functions respectively, are shown (curves 
with symbols (red)), with the dashed one from \tpTi)- The mean-field results that \ipTi) 
and \1pT2} yield are also shown (no symbols). We see that the two sets of CPMC results, 
although not in exact agreement, over-correct the mean-field results and show consis- 
tency. The statistical error bars on the CPMC results are indicated by small dots, the 
maximum being several times the symbol size. 




interaction strength is increased to 4. The statistics is enough to discern in 
this case the irregular oscillations of Dd(l) around zero. 

Calculation was also done for the 12 x 12 system at U = 8 with the 
same (free-electron) trial wave function. The error bars become substan- 
tially larger. Dd(l) oscillates around zero and shows no indication of en- 
hancement at large distances [2]. 

We also note that the behavior of Dd(l) in Fig. 4 at short distances 
is just the opposite of that at large distances, i.e., Dd(\) increases as U is 
increased. Due to the large discrepancy in magnitude between short and 
long distances, the integrated value of Dd(l) would therefore not give the 
correct trend even at system sizes as large as these. 

In summary, the near-neighbor d x 2_ y 2-wave electron pairing correlation 
seems to disappear as the interaction strength increases and no indication 
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0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 1 1 .0 

l/l 

Figure 4- d x 2_ y 2-wa,ve pairing correlation Dd(l) vs. distance I = |1|. TOP: Da(l) for a 
0.85-filled 12 x 12 lattice at interaction U — 2 and 4. The inset shows the the portion 
with I > 3. Shown together is the ground-state pairing correlation for the non-interacting 
system, whose wave function is also the trial wave function used in the CPMC calculations 
for both values of U. At large distances (inset) Dd(l) systematically decays as U is 
increased. Also note that at the irrelevant (short) distances the behavior is the opposite. 
BOTTOM: Same plot, for a 16 x 16 system. Same trend as that of 12 x 12 continues, but 
now the long-range correlation vanishes at U — 4. 
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of enhancement is found over the non-interacting system. Indeed, the be- 
havior of the pairing correlation appears similar across a fairly broad range 
of electron fillings, including the half-filled case, which exhibits pairing cor- 
relations of similar magnitude [2]. 

We emphasize again that, due to the approximate nature of the CPMC 
method, the possibility can not be completely ruled out that the results were 
biased such that the pairing correlation was suppressed. However, the many 
benchmark calculations, and particularly the variety of self-consistency 
checks, indicate strongly that the results are accurate and robust. 

4. Finite-temperature CPMC algorithm 

Extension of the CPMC algorithm to finite temperatures will allow us to 
remove the second set of question marks in Table 1. 

While the standard method of Blankenbecler, Scalapino, and Sugar 
(BSS)[4, 27] has been widely applied, the exponential scaling has hindered 
its ability to effectively study true phase transitions. The difficulty with 
developing a finite-temperature counterpart of CPMC lies in the implicit 
nature of the path-integral picture in the BSS formalism. As discussed in 
Table 1, paths in the grand-canonical formulation, contrary to ground-state 
projector QMC, do not originate or end at a single explicit point in Slater 
determinant space. Instead they involve many points, due to the trace over 
TV (the number of particles, not the number of sites). Indeed the points do 
not even have the same "dimension" [4] . A naive application of the current 
way of imposing the constrained path approximation, which would require 
separate constraining wave functions and conditions for paths with different 
end points, would therefore not be practical. 

We have developed a mechanism to incorporate the constraint into the 
analytic evaluation of the trace [3], which we are currently testing. In Table 

5, we show preliminary results on total energies per site, for a simple case 
of the 4x4 system. The chemical potential is adjusted so that the average 
filling fraction corresponds to 5 f and 5 J. electrons. As high temperatures, 
the new algorithm reproduces the results from AFQMC, while at a low 
temperature of T = 0.0625 {(3 = 16), it is in good agreement with the 
ground-state result from exact diagonalization, as well as that from the 
T = OK CPMC algorithm (cf Table I). 
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TABLE 5. Preliminary data from the new fi- 
nite-temperature algorithm. The computed total 
energies per site, as a function of temperature, are 
compared with results from AFQMC[34] or exact 
diagonalization. The system is 4 x 4 at U = 4. The 
chemical potential is adjusted so that the average 
filling fraction corresponds to 5 | and 5 J. electrons. 
Statistical errors are in the last digit and are in- 
dicated in parentheses. The * indicates the exact 
ground-state energy from exact diagonalization. 



T 


1 


0.25 


0.0625 


CPMC 


-0.8459(2) 


-1.1947(6) 


-1.2237(6) 


AFQMC 


-0.8457(1) 


-1.1947(2) 


-1.2238* 
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